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A deeper understanding of nonequilibrium phenomena is needed to reveal the principles govern- 
ing natural and synthetic molecular machines. Recent work has shown that when a thermodynamic 
system is driven from equilibrium then, in the linear response regime, the space of controllable 
parameters has a Riemannian geometry induced by a generalized friction tensor. We exploit this ge- 
ometric insight to construct closed-form expressions for minimal-dissipation protocols for a particle 
diffusing in a one dimensional harmonic potential, where the spring constant, inverse temperature, 
and trap location are adjusted simultaneously. These optimal protocols are geodesies on the Rie- 
mannian manifold, and reveal that this simple model has a surprisingly rich geometry. We test these 
optimal protocols via a numerical implementation of the Fokker-Planck equation and demonstrate 
that the friction tensor arises naturally from a first order expansion in temporal derivatives of the 
control parameters, without appealing directly to linear response theory. 

PACS numbers: 05.70.Ln, 02.40.-k,05.40.-a 



I. INTRODUCTION 

There has been considerable progress in the study of 
nonequilibrium processes in recent years. For example, 
fluctuation theorems relating the probability of an in- 
crease to that of a comparable decrease in entropy during 
a finite time period have been derived [IH5] and exper- 
imentally verified [5H5] in a variety of contexts. More- 
over, other new fundamental relationships between ther- 
modynamic quantities that remain valid even for sys- 
tems driven far from equilibrium, such as the Jarzynski 
equality [TDHT5] . have also been established. Interest- 
ingly, some of these ideas were independently developed 
in parallel within the machine learning community |14j . 
as ideas from nonequilibrium statistical mechanics are in- 
creasingly finding applications to learning and inference 
problems p~5l [16] . 

For macroscopic systems, the properties of optimal 
driving processes have been investigated using thermo- 
dynamic length, a natural measure of the distance be- 
tween pairs of equilibrium thermodynamics states 
I22j . with extensions to microscopic systems involving a 



metric of Fisher information [331 IMj ■ Recently, a linear- 
response framework has been proposed for protocols that 
minimize the dissipation during nonequilibrium pertur- 
bations of microscopic systems. In the resulting geomet- 
ric formulation, a generalized friction tensor induces a 
Riemannian manifold structure on the space of param- 
eters, and optimal protocols trace out geodesies of this 
friction tensor 25j. 

In this article, we make use of Riemannian geometry 
theorems to simplify the problem of optimizing proto- 
cols. To illustrate the power of these geometric ideas, we 
consider a simple, but previously unsolved, stochastic 
system and calculate closed-form expressions for optimal 
protocols. We test the accuracy of our approximation 
by numerically comparing our optimal protocols against 
naive protocols using the Fokker-Planck equation. We 
conclude by demonstrating that our inverse diffusion 
tensor framework arises naturally from a first order 
expansion in temporal derivatives of the control parame- 
ters, without appealing directly to linear response theory. 



II. DERIVATION OF THE EXCESS POWER 
FOR VARIABLE TEMPERATURE 



* pzulkowski@berkeley.edu 

t dasivak@lbl.govl Present address: Center for Systems and Syn- 
thetic Biology, University of California, San Francisco, CA 94158; 
david.sivak@ucsf.edu 

t gcrooks@lbl.gov 

§ deweese@berkeley.edu 



For a physical system at equilibrium in contact with 
a thermal bath, the probability distribution over mi- 
crostates x is given by the canonical ensemble 



tt(x|A) = exp/3[F(A) - E(x,X)] 



(1) 
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where j3 = (/cbT) -1 is the inverse temperature in natural 
units, F(X) is the free energy, and E(x,X) is the system 
energy as a function of the microstate x and a collection 
of experimentally controllable parameters A. 

In equilibrium, the thermodynamic state of the system 
(the probability distribution over microstates) is com- 
pletely specified by values of the control parameters, but 
out of equilibrium the system's probability distribution 
over microstates fundamentally depends on the history of 
the control parameters A, which we denote by the con- 
trol parameter protocol A. We assume the protocol to 
be sufficiently smooth to be twice-differentiable. 

The usual expressions for heat and work [2l)H2"9"] as- 
sume that the temperature of the heat bath is held con- 
stant over the course of the nonequilibrium protocol. 
Following the development of methods to handle time- 
varying temperature described in section 1.5 of [30], and 
preceding Eq. (4) of [31], we argue that the unitless en- 
ergy f3E{x, A) (normalized by the natural scale of equilib- 
rium thermal fluctuations, k^T = /J -1 , set by equiparti- 
tion) is the fundamental thermodynamic quantity. Thus 
when generalizing to a variable heat bath temperature, 
we arrive at the following definition for the average in- 
stantaneous rate of (unitless) energy flow into the system: 



dt 



f3E(x,\) 



(2) 



where angled brackets with subscript indicate a nonequi- 
librium average dependent on the protocol A. For con- 
stant /3, this reduces to the standard thermodynamic def- 
inition |25j . With this definition, we can prove that for 
systems obeying Fokker-Planck dynamics, excess work 
is guaranteed to be non-negative for any path, which is 
not true of the naive definition (see § |III[ ). Nonetheless, 
a deeper understanding of the subtleties involved in our 
modified energy flow definition (Eq. ([2])) calls out for fur- 
ther study. 

Eq. ([2]) may be written as 



/ dx T d (J3E) 
\~dt 



dx 



(x,X) 



d\ T d (J3E) 



dt 



ox 



Or, A) 



The first term represents energy flux due to fluctuations 
of the system at constant parameter values and naturally 
defines heat flux for nonequilibrium systems. The second 
term, associated with an energy flux due to changes of 
the external parameters, defines nonequilibrium average 
power in the general setting of time-variable bath tem- 
perature. 

The average excess power exerted by the external agent 
on the system, over and above the average power on a 
system at equilibrium, is 



Here X = — d< "g^ are the forces conjugate to the control 



dX 1 

~dt 



(AX) 



(4) 



parameters A, and A_X"(i ) = X(t ) — (X) X (t a ) * s the 
deviation of X(to) from its current equilibrium value. 
Applying linear response theory 32 , 



(AX(t )) A 



X (t - t') ■ [A(to) - A(t')] dt' (5) 



where Xij(t) = — dll^j(t)/dt represents the response of 
conjugate force Xi at time t to a perturbation in control 
parameter \ 3 at time zero, and 

^\t) = {8X i {Q)SX i (t)) x{!tor (6) 



For protocols that vary sufficiently slowly [35] , the result- 
ing mean excess power is 



dX 1 

~dt 



fl(A(*o)) 



dX 

~dt 



(7) 



for inverse diffusion tensor 
<?y = /?(to)Cy(A(io)) = 



dt'^Xj (0)SXi(t')) 



Mto) 



(8) 

where Qj is the friction tensor in control parameter 
space from [2SJ. We will construct geodesies using this 
inverse diffusion tensor gij. 



III. THE MODEL SYSTEM AND ITS INVERSE 
DIFFUSION TENSOR 

We consider a particle (initially at equilibrium) in a 
one-dimensional harmonic potential diffusing under iner- 
tial Langevin dynamics, with equation of motion 

my + k(t) (y - y {t)) + ( c y = F(t) , (9) 

for Gaussian white noise F(t) satisfying 

<F(t)> =0, (F(t)F(t>)) = ^ ) 5(t-t>) . (10) 

Here £° is the Cartesian friction coefficient. We take as 
our three control parameters: the inverse temperature 
of the bath /3, the location of the harmonic potential 
minimum yo, and the stiffness of the trap k [see Fig.Qa)]. 
The conjugate forces are 



13k (y - y Q ) , • 



2m 



^{y-yo) ,-^{y-yo) 



(ii) 

This model can be experimentally realized as, for in- 
stance, a driven torsion pendulum [33l [34], 
The excess work 



((Wcx) 



dt P(t)V m (t) 



(12) 
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(a) 



Large k High /3 



Small k Low /3 




FIG. 1. (a) Our model system. A particle (black dot) diffus- 
ing in a harmonic potential with adjustable spring constant k, 
position i/o, and inverse temperature /3 = (indicated by 
thermometers), (b) Representative optimal protocols (orange 
and red curves) plotted for two of the three control parame- 
ters, k and f3. An optimal protocol (e.g., red curve) results 
in the minimum dissipation for any path taking the system 
from one particular state (black square) to another (black tri- 
angle) in a fixed amou nt o f time, (c) A change of variables 
{P,k} — > {z, x} (Eq. ( 34 1 ) reveals that our model system 
has an underlying structure described by hyperbolic geom- 
etry, represented here as the Poincare half-plane, in which 
geodesies form half circles (orange curve) or vertical lines (red 
line), (d) A piece of the (z,x) manifold may be isometrically 
embedded as a saddle in M 3 . The distortions in each of these 
two optimal paths as shown in panels (b) and (c) reflect the 
curvature of this manifold. 



is non-negative. Assuming the system begins in equilib- 
rium, the relative entropy D\f || 7r(x|A(i))] corresponds 
to the available energy in the system due to being out of 
equilibrium |35j , and bounds the excess work from below. 
Here, tt is the equilibrium distribution (Eq. ([IJ) defined 
by parameters X(t), and / = f(y,p,t) is the nonequilib- 
rium probability distribution. The time derivative of the 



relative entropy may be written as 



;D[f\\n(x\\(t))] 



Of 



dt lJ " " v '" W ' J ./ dt 
which follows from the identity 



log(f/n(x\X(t)))+P(t)V ex (t), 
(13) 



/ 



dX 1 

~dt 



(AX(t)) A . (14) 



The first term of Eq. (131 simplifies to 



C c / -en 
J 1 ' 



1 



d 



f \dp 



e2m / 



< o. 



(15) 



Integrating Eq. ( 13 1 from to to proves the relative en- 



tropy bounds the excess work from below. Since this 
quantity is always non-negative, so is the excess work; in 
fact, for any finite-duration path visiting more than one 
point in parameter space, it is strictly positive, yield- 
ing a well-behaved metric in our geometrical formalism. 
See [35] for related results in the special case of constant 
temperature. Note that, unlike our modified definition 
for work, the naive definition /J dtV ex {t) may be nega- 
tive for particular protocols that vary /?. 

Calculation of the time correlation functions in Eq. ^ 
requires knowledge of the dynamics for fixed control pa- 
rameters. We may write any solution to the equation 
of motion as a sum y^ + y p of a homogeneous part yh, 
which depends on the initial conditions and is indepen- 
dent of F(t), and a particular part y p , which has vanish- 
ing initial conditions but depends linearly on F(t) (see, 
for instance, Theorem 3.7.1 in |37jV Explicitly, we may 
write 



Vp(t) 



( V l H\s)vF\t)-y£\t)yj?\s) \ Fjs) 



ds 



(16) 

where y^\t) for i — 1, 2 are independent solutions of the 
homogeneous equation. It follows immediately that 



(2) 



(17) 



where the constants C\ , C2 are determined by initial con- 
ditions. 

For Gaussian white noise F(t), it is easy to show that 
the particular piece y p does not contribute to the equi- 
librium time correlation function (SXj(Q)SXi(t)) . For 
simplicity and without loss of generality, consider the 
correlation function (^y(t) 2 <5?/(0) 2 ). Expanding this ex- 
pression, 

(Sy(tf6y(0) 2 ) = (y(tfy(0f)-(y(tf)(y(0f) , (18) 
and substituting y(t) = yh(t) + y p (t), we find 

(Sy(t) 2 5y(0) 2 ) = (y h (t) 2 y(0) 2 ) + (y P (tfy(0) 2 ) 

-(vh(tf)(y(0f)-(y P (tf)(y(0) 2 ) (19) 

+ 2( (y h (t)y p (t)y(0) 2 ) - (y h (t)y p (t)) (y(0) 2 ) ) . 
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Angled brackets denote an average over noise and initial 
conditions. 



Integrating these expressions, we obtain 



According to Eq. (16), the particular solution y p does 
not depend on the initial conditions. It follows immedi- 
ately that 

(y P (t) 2 y (o) 2 ) - (y P (t) 2 ) (y(0) 2 ) = o . (20) 

Furthermore, since y^ depends only on the initial condi- 
tions and is independent of the noise, 

(y h (t)y P (t)y(0) 2 ) - (y h (t)y p (t)) (y(0) 2 ) = , (21) 

follows from the assumption that (F(t)) = 0. To sum- 
marize, 



(Sy(t) 2 Sy(0) 2 ) = (Sy h (t) 2 Sy(0) 2 ) . 



(22) 



For each of the time correlation functions needed to com- 
pute the inverse diffusion tensor, it is generally true that 
yh(t) may be substituted in the average for y(t). 

Without loss of generality, let us assume for the mo- 
ment that (C c ) - 4fcm > 0. If we define 



2m 2 V V m / m 



(23) 



then the homogeneous solution with initial conditions 
{?/(0),p(0) = my(0)} is given by 



Vh(t) = y 



p(0) + mr_ (y(0) - y ) ^_ r+t 
m(r_ — r+) 
p(0) + mr + (y(0) - y ) _ t 
m(r + — r_) 



(24) 



where yo is the fixed trap position. For convenience, let 
us define Y = y~yo- Assuming that the initial conditions 
{y(0),p(0)} are distributed according to the equilibrium 
Boltzmann distribution 7r[y(0),p(0)] oc e -PE[y(o), P (o)] for 

E[y,p] — 2m + \kY 2 , we obtain the following identities: 



m®sY'{o)) = — -j- - 2 

(«P) (r + -r-) 



(25a) 



(SY 2 (t)SY 2 (0)) 
(SY 2 (t)6p 2 (0)) 
(SY 2 (t)Sp 2 (0)) = 



2rir 2 



-rit „— r_t\2 



fc 2 /3 2 (r + — r_) 
2 

/3 2 (r + — r_)^ 
2 

/3 2 (r+-r_) 2 



e 'f — e 



e '-t- — e 



(25b) 



(25c) 

^ — T-t „ „ — r+t \ ^ 



(25d) 



dt(SY 2 (t)6Y 2 (0)) 



dt(5Y 2 (t)SY 2 (0)) = 



dt(SY 2 (t)Sp 2 (0)) = 
dt(5Y 2 (t)S P 2 (0)) = 



(C c 



fc 2 /3 2 C c 
1 



fcm 



m 

27c 



t 2 

m 

C^ 2 



Thus the inverse diffusion tensor is 
/^/3 

?Ti. " 



4C C 



V 



/3 2 l * ^ fern J /3fc \ A ^ km J 



J_ (0+ KU- , - 

/9fc \ ~ fcm J fc 2 



/cm 



(26a) 
(26b) 
(26c) 
(26d) 

\ 

/ 

(27) 

which endows the space — oo < yo < oo, < /3 < oo, < 
fc < oo with a Riemannian structure. 



IV. BRIEF REVIEW OF RIEMANNIAN 
GEOMETRY 

We recall some definitions from Riemannian geometry 
and establish notation (see [38H40] for details). For a 
smooth Riemannian manifold M endowed with metric 
tensor g, the Christoffcl symbols arc defined as 

Y\ k = Ig* (d i9kl + d k g u - d l9ik ) (28) 

where g % i denotes the matrix inverse of the metric. We 
employ the Einstein summation convention here (and as- 
sume it throughout). The Riemann tensor, constructed 
from the Christoffel symbols, measures the curvature of 
the manifold M and is given in local coordinates by 

R l jM = d k ^ l ji - dir l jk + r^r l mfc - r™ k r l ml . (29) 

Contracting indices gives the Ricci tensor iiy and the 
Ricci scalar R, 



Rij — R uj , R = g lJ R 



13 > 



(30) 



which are useful for determining the curvature content of 
the manifold M . Geodesies are defined in local coordi- 
nates by 



d 2 X i , d\ j d\ k n 

-nr + r ik-, — -r- = • 
dr 2 J dr dr 



OPTIMAL PROTOCOLS 



(31) 



Though one can write down the geodesic equations for 
the metric Eq. (27) in the (yo,/3,k) coordinate system, 
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more insight is gained by finding a suitable change of 
coordinates. Consider the lower right 2x2 block of the 
metric Eq. ( |27[ ) which is the metric tensor for the two- 
dimensional (p, k) submanifold. A direct calculation of 
this submanifold's Ricci scalar yields R = —2( c /m which 
is constant and always strictly negative. 

Theorems from Riemannian geometry |38j imply that 
this constant negative-curvature submanifold is isomet- 
rically related to the hyperbolic plane. In our construc- 
tion, we choose the Poincare half-plane representation 
of the hyperbolic plane, which is described by {(z, x) £ 
M. 2 ,x > 0} with metric tensor given by the line element 

ds 2 = g i jdx l dx : > — dx + 2 dz . The geodesies of the hyper- 
bolic plane (see Fig. [l| are half-circles with centers on 
the z-axis and lines perpendicular to the z-axis. Fig.JIJc) 
shows two geodesies in (z, x) coordinates. The portion of 
the hyperbolic plane {(z,x) £ 
be isometrically embedded in 



2 ,x > l,z £ [0, 7r]} may 
' using the map 



— cos z, log 



(V* 




(32) 



The geodesies of Fig. [T] (c) and the part of the hyperbolic 
plane containing them are embedded in K 3 in Pig. 13d). 

The line element associated with the submanifold met- 
ric tensor, 



ds 2 = 




(33) 



km 



dk 1 



is coordinate-invariant since it measures geometric dis- 
tances. Thus we may construct an explicit coordinate 
transformation, 



1 



1 



(34) 



to demonstrate the equivalence of the submanifold with 
a portion of the Poincare plane. Note that x is propor- 
tional to the classical partition function of the system 
in equilibrium, and z is proportional to the equilibrium 
variance of y — y . Inverting Eq. (34 1, and substituting 



into Eq. ( 33 ) gives the metric tensor in (z, x)-coordinates, 
m dx 2 



ds' 



dz 2 



(35) 



The line element corresponding to the metric of the 
full three-dimensional manifold in Eq. (27) is 



2 _ mz dy 2 



ds 2 = 



dx 2 + dz 2 



(36) 



in (y ,z,x) coordinates. To fully exploit the machinery 
of Riemannian geometry to find closed-form geodesies, 



we look for Killing fields of Eq. ( 36 1 . In general [551 - 



[40] . isometries of a metric are generated by the Killing 
vector fields K which are themselves characterized by the 
Killing equation 



2Y%K k = . (37) 



While directly solving this system of equations may be 
difficult, certain characterizations of Killing vectors help 
circumvent this difficulty. For instance, if in a given co- 
ordinate system the metric tensor components are inde- 
pendent of a coordinate x l , then the coordinate vector 
d x % is a Killing field [ID]. Hence, d Vo is clearly a Killing 
vector field. Examining the full set of Killing equations 
shows that 



K = y d yo + 2xd x + 2zd z 



(38) 



is also a Killing vector field. There may be more solutions 
to the Killing equation yet to be discovered. 

In general [40 , for Killing vector Ki the quantity K t 



is conserved along the geodesic described by A. 
follows from 



d 



K, 



d\ l 
dr 7 



ViK 3 



d\ l dX 3 
dr dr 



1 dr dr 



This 



(39) 



The first term of the equation vanishes by the definition 
of the Killing field and the second term vanishes by the 
geodesic equation Eq. (31 1. For the three-dimensional 



inverse diffusion tensor, we have the following two con- 
served quantities associated with Killing fields: 

z(t) dy 2 dx z(t) ,dy 2z(t) dz 

+ -^vo(t)-t- + -tt-.-t- (40) 



x 2 (t) dr ' x(t) dr x 2 (t) 



dr x 2 (t) dr 



To solve the geodesic equations, note that the velocity 
of the geodesic (i.e., its tangent vector) must have con- 
stant norm [3"5H4"U] . For convenience, we choose the norm 
so that 



c\ x 2 (t) 




(41) 



where we have used the first conserved quantity of 
Eq. ( 40 1 . We combine this with the full geodesic equation 
for x(t), to decouple x(t) from yo(r) and z[t): 



d 2 x 2 
dr 2 x(t) 

which has solution 



dx \ ' 
dr) 



+ x(t) = 



x(t) = r sech(r) 



(42) 



(43) 



When z(t) is constant, the geodesic equation for z implies 
that ?/o is also constant, giving a geodesic straight line in 
the constant-z submanifold. 



When z(t) is not constant, Eqs. (41 1 and (43) imply 

r~ r '■{ t i 

(44) 



x\t) 



dz\^ | c{ x 4 (t) 



dr J 



z(r) 
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which integrates to 

z(r) = hT x [c2 — rtanh(r)), 

where 



(45) 



h(£) = £Jl- f+|c?log 2£(l+Wl-f)-c? • (46) 



The Killing conserved quantities of Eq. (|40j), together 
with x(t) and z(t), yield 



2/o if) = -B-ci log 



cf + 2/i 1 (c 2 - rtanh(r)) x 



/i 1 (c 2 — r tanh(r)) 



(47) 



Let (yo,i, Xj, Zi) and (yo t f,Xf, Zf) denote the endpoints 
of the geodesic. Define AA = A/ — Aj and A = A, + A/ 
for A e {yo, z}. Defining h = !^Zi+Mgtj ancl A/i = 
/i (zy) — ft, (zj) , the constant c 2 may be written as 

- _Ax 
c 2 = h + x-— 
Ah 

and r is given by 
The constant is given by 



(48) 



(49) 



E = y ,i + ci log ( -c? + 2zi | 1 + \j 1 - J 



and ci is determined by the equation 




(50) 



Ay = -ex log 



2z, 1 



log -c? + 2zAl 



The parameter r ranges between the values 

Ax 



and 



/ , Al _ \ i fXi\ 

Ti = sgn I An + * l ~£j l x ) secn ^ — J 
r / = -sgn(Afc-2^x) sech" 1 (^) . 



(51) 



(52) 



(53) 



When y is held fixed, the geodesies are precisely 
those of the hyperbolic plane as expected. Furthermore, 
these geodesies are necessarily minimizing by virtue 
of the constant, negative Ricci scalar [38j [39]. Several 
example geodesies are displayed in Fig. [2] 



VI. COMPUTING DISSIPATION 
NUMERICALLY 



We validate the optimality of these geodesies by calcu- 
lating excess work directly from the Fokker-Planck equa- 
tion. In full generality, the mean excess work as a func- 
tional of the protocol X(t) = (yo[t), j3(t),k(t)) is 



dt f3V ex (54a) 
^ dt (f 2^ + l {iy ~ Vo)2) + kP ) 



+ fc/3yo(yo - y) 



(3 2k 



(54b) 



Here angled brackets denote averages over the nonequi- 
librium probability density f{y 1 p,t). 

Standard arguments [32] yield the Fokker-Planck equa- 
tion for the time evolution of f(y,p,t), 



df Vdj df C c d[pf] C d 2 f _ 

dt + mdy~~ m [V ~ Mt)l fy~m^p~~W)W ~ ' 

(55) 

By integrating Eq. ( 55 ) against y, p, etc., we find a system 




10 

8 '-_ 
6 1 
4 1 
2 1 




FIG. 2. Optimal protocols differ substantially from linear in- 
terpolation (red dashed lines). Blue solid curves represent 
geodesies of the inverse diffusion tensor, and are thus opti- 
mal protocols for transitioning the system from one state to 
another in a fixed amount of time. Blue dots indicate points 
separated by equal times along each of the eight optimal paths 
shown. 
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of equations for relevant nonequilibrium averages: 



(y) 
(p) 



dt 
d 

dt 
Jt {py) 



(p) 



m 



(p 



m 

2\ 



(p) -k(y- yo) 



m 
2 



k(y 2 



C 



— (py) + k vo(y) 

to 



^(y ) = —(py) 



dt 



to p 



(56a) 
(56b) 
(56c) 
(56d) 
(56e) 



Following the derivation of the friction tensor in [25] 
would require us to use linear response theory and to 



supplement the system Eq. ( 56 ) by initial conditions 



(y) (0) 


= yo(Q) 


(57a) 


(p) (0) 


= 


(57b) 


(y 2 ) (0) 


= yo(0)2+ fc(OMO) 


(57c) 


(py) (0) 


= 


(57d) 


(p 2 ) (0) 


TO 


(57e) 



We solve these equations numerically and compare a 
geodesic protocol with naive protocols in Fig. [3] 

This system has three natural dimensionlcss quantities 



C c Ai kAt Z 2 m 2 /? 



(58) 



dependent upon characteristic scales for (inverse) tem- 
perature /3, length I, spring constant k and the protocol 
duration At. These suggest at least two plausible mea- 
sures of distance from equilibrium [25] . A corresponds to 
the ratio of two timescales, the timescale p for frictional 
damping and the timescale of the perturbation protocol 
At. Likewise, B is the ratio of two powers during changes 
of yo, the dissipative power ( c (Ayo/At) 2 and the elastic 
power k(Ay ) 2 /At. As A decreases and as B decreases, 
the system will remain closer to equilibrium during the 
course of the nonequilibrium perturbation, and hence our 
near-equilibrium approximation will be more accurate. 

This intuition is confirmed in our numerical calcula- 
tions: with A <C 1 and B <C 1, the dissipation of geodesic 
protocols obtained numerically via Fokker-Planck agrees 
with the inverse diffusion tensor approximation to better 
than .1% (see Fig. [3]). Note that, while the inverse 
diffusion tensor approximation is excellent for optimal 
protocols and small deviations thereof, it can deviate 
substantially from the exact result for large deviations 
from the geodesic. 



VII. THE INVERSE DIFFUSION TENSOR 
ARISES NATURALLY FROM THE 
FOKKER-PLANCK EQUATION 

If we neglect terms involving derivatives of protocols 
of degree two and higher, we may find an approximate 
solution to the Fokker-Planck system: 
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Substituting this into the expression for mean excess 
power Eq. (54b I, we recover Eq. 0. The argument 
above suggests that the emergence of the inverse diffu- 
sion tensor from the Fokker-Planck equation may follow 
from a perturbation expansion in small parameters. 



VIII. DISCUSSION 

We have employed geometric techniques to find op- 
timal protocols for a simple, but previously unsolved, 
stochastic system. Calculation of the Ricci scalar for 
a submanifold pointed to a change of coordinates that 
identified the submanifold with the hyperbolic plane 
and greatly simplified the metric for the full three- 
dimensional manifold. This simplification, combined 
with the identification of a Killing field, permitted cal- 
culation of an exact closed-form expression for geodesies. 
Exact calculations using the Fokker-Planck equation con- 
firmed that geodesies in the (/3, k) -submanifold do indeed 
produce less dissipation than any comparison protocol we 
tested. 

In addition to being useful for identifying optimal pro- 
tocols, we expect that the Ricci scalar will turn out to 
have an important physical interpretation. Riemannian 
geometry has been useful for the study of thermodynamic 
length of macroscopic systems jHJ 02], and there has 
been some speculation about the role of the Ricci scalar in 
that setting I42| , but the interpretation of R arising from 
the inverse diffusion tensor remains ambiguous. We hope 
that further study of these geometrical ideas extended to 
nonequilibrium systems will help clarify its role. 

It would also be interesting to establish a physical 
interpretation for the conserved quantities arising from 
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FIG. 3. Geodesies describe protocols that outperform naive 
straight line paths in parameter space. A geodesic between 
two fixed points in the (/3, fc)-plane (black) and several com- 
parison protocols are pictured above. The comparison pro- 
tocols were generated via a linear interpolation between the 
constant speed straight line (pink) and the geodesic. The tick 
marks represent points separated by equal times. The solid 
pink dots correspond to the constant speed parametrization 
of the line whereas the open red circles correspond to the 
optimal parametrization along this straight path. The ra- 
tio of excess work to that of the geodesic protocol is: 6.12 
(pink circle), 4.37 (red open circle), 3.67 (cyan downward- 
pointing triangle), 1.38 (orange upward-pointing triangle), 
1.00 (black star (geodesic)), 2.86 (magenta square), and 6.29 
(green circle). These ratios are plotted in the inset figure 
along with a graph of the ratio as a function of the interpolat- 
ing parameter (light gray curve). All excess work values were 
calculated using the Fokker-Planck system Eq. (56 1. Here, 
A = 10~ 2 ,B = 10~ 3 ,A/ = 1, placing the system within the 
near-equilibrium regime and ensuring accuracy of the inverse 
diffusion tensor approximation. 



Killing fields in this context. We found two conserved 
quantities (see Eq. (40)), which may be the only ones, 
but this model could have as many as six, given that 



there might be as many as six unique globally smooth 
Killing fields for this three-dimensional model system. 
(In general, there are at most |n(n + l) independent 
globally smooth Killing fields where n is the dimension 
of the manifold [40].) 

In the course of developing our framework, we encoun- 
tered four distinct measures of the departure from equi- 
librium. The first two were dimensionless parameters, 
A and B, which have relatively straightforward physical 
interpretations — the timescale for frictional dissipation 
relative to the protocol duration and the ratio of the dissi- 
pative power to elastic power, respectively (see discussion 



The third was the disagreement between dissipation 
computed assuming linear response theory and the true 
dissipation. Empirically, we found that our linear re- 
sponse approximation was consistently accurate for all 
parameter regimes we tested in which both dimensionless 
parameters A and B were small, at least for protocols not 
too far from geodesies. Conversely, the linear response 
approximation appeared to break down for many cases 
we tested with at least one of these parameters of or- 
der unity or greater. However, the full extent of validity 
of the linear response approximation is not clear to us, 
suggesting an important direction for future research. 

Finally, we found that truncating to first order in tem- 
poral derivatives of the control parameters in our model 
was sufficient to yield the same inverse diffusion tensor 
formalism we originally derived using linear response the- 
ory. While it is plausible that these two types of linear 
approximations are directly related, further exploration 
is needed to uncover the relationship between linear re- 
sponse theory and truncating the model equations to first 
order in temporal derivatives. 

Our results are novel in three distinct ways. First, we 
included /3 as a control parameter, which is a natural 
extension of thermodynamic length {e.g. [4TJ [43] ) that 
is amenable to direct experimental confirmation. Our 
work generalizes the construction of |25j and opens up 
new experimental avenues for testing the validity of the 
framework. 

Secondly, our geodesic protocols optimize dissipation 
for simultaneous variation of all three adjustable param- 
eters; to our knowledge, no previous study has reported 
optimal protocols for any model system with three con- 
trol parameters. In [331 [35], Seifert and coworkers ele- 
gantly derived the exact optimal protocols for perturb- 
ing the position tjq and spring constant k separately, for 
both over-damped and under-damped Langevin dynam- 
ics. In [35] , Aurell and coworkers discuss the simulta- 
neous variation of the stiffness and the location of the 
trap. We note that our method misses the protocol jumps 
found in their analysis due to our smoothness assump- 
tions on the protocols. When this restriction on the 
differentiability of the curve is imposed, we found that 
any component of the optimal protocol (yo(t), P(t), 
generically depends on all components of both endpoints 
due to the non-trivial geometry of the parameter space. 

Finally, we successfully brought the machinery of 
Riemannian geometry to bear on a small-scale, nonequi- 
librium thermodynamic problem, revealing a surprisingly 
rich geometric structure. Concepts such as Killing vector 
fields, coordinate invariance and the Ricci scalar proved 
indispensable in the construction of optimal protocols. 
These results are encouraging and this approach may 
prove useful for understanding the constraints on the 
non-equilibrium thermodynamic efficiency of biological 
and synthetic molecular machines. 
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